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Abstract 



We propose a general formulation of simplicial lattice gauge theory inspired by 
the finite element method. Numerical tests of convergence towards continuum 
I— I results are performed for several SU{2) gauge fields. Additionaly, we perform 

"t^ simplicial Monte Carlo quantum gauge field simulations involving measurements 

'—H of the action as well as differently sized Wilson loops as functions of j3. 
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J> 1. Introduction 

in 

f^ General introduction. Gauge quantum field theory (QFT) has been extremely 

^^ successful in modeling the behaviour of fundamental high energy particle physics. 

" [ This is done using the standard model of particle physics, which is based on the 

C — ■ gauge symmetry group Q = C/(l)©5'f7(2)®5'?7(3). Quantum gauge field theories 

^^ based on such noncommutative gauge groups are also called Yang-Mills theories 

^, [B El El H] . Despite the massive successes of this model, there are still large 

. . difficulties in calculating low energy properties of quarks and gluons. When 

restricting to these quantum fields, the standard model reduces to the theory of 
Quantum Chromodynamics (QCD), with gauge group SU^Z). The problems is 
%2^ that through the effect of renormalization, the QCD coupling constant increases 

^ as interaction energies is decrease, in such a way that perturbation theory breaks 

down. This phenomenon is the source of confinement in QCD. Direct paper- 
and-pen calculation of masses and interactions among low energy bound states 
of quarks is therefore quite problematic. 



Lattice gauge theory. By discretizing QCD onto a lattice, a lot of these diffi- 
culties are removed. Lattice gauge theory (LGT) [51 [6] has proven itself to be 
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a powerful method of doing nonperturbative gauge theory calculations. It has 
therefore been, still is, and will for a long time be immensely useful in testing 
QCD against experimental results at low energy. 

Usually LGT models are formulated using a hypercubic lattice on a euclidean 
spacetime. Such a mesh preserves some discrete subgroups of the translational, 
mirror and 4d rotational symmetries. Note that a clever way of retaining con- 
tinuous symmetries while working on a lattice is to use random lattices [3 |Hl E] ■ 

The models are almost always defined so as to also preserve a discrete gauge 
symmetry. This has the beneficial effect of enforcing a vanishing gluon mass in 
the discrete model. 

Simplicial lattices. Simplicial meshes have been used for QCD simulations be- 
fore [ini El HH [131 im im mi, with promising numerical results. Here, we 
construct a simplicial gauge theory (SGT) based on the general mathematical 
concept of a simplicial complex, while preserving gauge invariance. This allows 
us to define SGT on a very general class of meshes, without restricting ourselves 
to a particular type of simplicial lattice. 

The construction of the gauge invariant SGT action functional is inspired 
by the finite element method (FEM) most commonly used for solving partial 
differential equations, particularly on complicated domains pTlfTBl [TOll^ CT| . 
The formalism therefore includes the use of finite element function spaces on 
simplicial meshes, and the concept of mass matrices. The latter has nothing to 
do with physical particle masses, and is therefore not to be confused with the 
usual mass matrix of quantum states within QFT. 

Through the use of the FEM formulation, and the massive resources of 
methods available within that subject area, we hope to gain advantages for 
QCD simulations in future implementations, in particular with regards to the 
possibilities of grid refinement. This could be useful in modeling some QCD 
phenomena, e.g. for highly concentrated gluon flux tubes between quarks where 
an increase lattice resolution might be desired. 

Computer simulation. The mathematical proof of consistency between the SGT 
and continuous Yang-Mills gauge theory action is described in a companion 
paper [22] , along with a description of the more comprehensive Yang-Mills-Higgs 
model. In the current article we are content to provide numerical evidence for 
convergence towards exact continuum results for several choices of gauge field 
configurations. In addition, we perform Monte Carlo quantum pure gauge field 
theory simulations for the gauge group SU{2) in temporal gauge, as a proof-of- 
concept for SGT. Observable measurements include expectation values of the 
action density as well as a series of different Wilson loops. 

Outline. Section [2] contains a short repetition of the fundamental definitions of 
gauge symmetry and the continuous spacetime Yang-Mills action in subsection 
|2.1[ the basics of traditional lattice gauge theory in subsection |2.2[ as well as 



an introduction to the proposed SGT action in subsection 2.3 In section 3.1 



we report on the numerical convergence of the SGT action towards the exact 



continuum value for several different cases of SU{2) gauge fields, as well as 
similar results from traditional LGT. Theoretical results proving convergence 



for general gauge fields can be found in [2^. In section 3.2 we perform Monte 
Carlo quantum field theory simulations in order to observe that SGT correctly 
reproduces the basic aspects of the SU{2) quantum field theory. We draw 
our conclusions in section |4j [Appendix A| contains a short introduction to 
elementary aspects of simplicial complexes, and some notes about basis functions 



and mass matrices that are use in our construction of SGT.[A"ppendix^ contains 



a calculation of strong and weak coupling limits for a Wilson triangle and the 
action density. Lastly, [Appendix C| contains a short discussion of some aspects 
of the numerical computer implementation. 

2. Construction 

2.1. Continuous gauge theory 

Consider the spacetime domain M = M x S*, where M is time and S C 
K^. The domain M represents cither lorentzian or euclidean spacetime, in each 
case equipped with the appropriate metric. In the standard orthonormal M- 
basis {ep}^=o,i,2,3, a general point a; G M has components {x'^}^=oa,2,3- Greek 
indices run from to 3, and Latin indices from 1 to 3. 

Furthermore, in this article we shall consider pure SU{2) gauge theory. How- 
ever, the construction presented is applicable to any gauge theory based on a 
compact Lie group Q which can be represented by a subgroup of the complex 
unitary n x n matrices. We define the real-valued scalar product on Q as 

g'-g:^mr{g'g"), (1) 

where g^ is the hermitian conjugate of a matrix g. 

The connection between the continuous theory and the discrete simplicial 
theory is most easily seen in a coordinate free formulation. Thus, we start with 
a coordinate free formulation, before we give the more familiar coordinate based 
one. 

The free variable in pure Yang-Mills theory with gauge Lie group Q is a gauge 
potential or more formally a one form A on M, with values in the corresponding 
gauge Lie algebra g. For simplicity of notation, we hereby specify Q = SU{2) 
and g = su(2). We split A into temporal and spatial components A = {Aq,A). 
In this context, Aq can be thought of as a scalar functioij^ , and A as a spatial 
vector. The curvature (field strength) of such a one form is given by 

F{A) - d^ + ^ [A, A] = doA + dAa + dA + i[A, A^^] + '- [A, A], (2) 

where d — {dQ,d), do and d denote exterior derivative in the temporal and 
spatial directions respectively, and [•,•] is the commutator between Lie algebra 

^However, not a scalar in the sense of spacetime symmetry transformation properties. 



valued one forms. We choose the basis {i°}o=i,2,3, where t°- :— cr°/2, for 5u(2), 
where {(J°'}a=i,2,3 are the Pauh matrices. Thus, we can expand the gauge field 
into components, A = A°'t°' . We also have 

[A, A] = ^ ^" A yl^[r , t^] = Yl i^^^^A" A A^t\ (3) 

ab abc 

where e'^^'^ is the antisymmetric Levi-Civita symbol with e^^^ = 1 and A is the 
wedge product (exterior product). For later convenience we split the curvature 
in a temporal and spatial part 

F\A)^doA + dAo + t[A,Ao], F%A)^dA+'-[A,A]. (4) 

The action that defines the gauge theory is the functional 

S[A] = -^ / \F{A)\^ = -^ / \F\A)\^ + \F%A)\\ (5) 

*6 JM *S JM 

where the norms are generated the metric and e is the dimensionless Yang-Mills 
coupling constant. 

A gauge transformation is defined by a choice of G{x) E SU{2) for each 
x g M, and transforms the gauge field as 

Aq >^ G {Aq + dt) G-\ A>->G{A + d)G-^. (6) 

Note that the action S[A] is invariant under such gauge transformations. For a 
more precise mathematical exposition, see |22| . 

A formulation more familiar within physics is obtained by expressing the 
one form and curvature in coordinates. In other words, one decomposes the 
one-form A"" in the basis {dx^}, i.e. A" = X^u^u'^^'^- The exterior derivative 
of such a one-form is given by 

dA" = J2 d^A^dx" Adx'' = J2 lidf^K - 5^^;)rfa;^ A dx" . (7) 

Furthermore, the curvature is given by F^ = ^ ^F^^dx^ A dx^ ^ where 

K^ = d,Al - d,A; - e'^'^A^Al. (8) 

Finally, the action can be expressed as 

S=^ fT.K^F''"'dx, (9) 

*" JM 

the usual coordinate dependent expression for the Yang-Mills action functional. 



2.2. Lattice gauge theory 

To see the connection between lattice gauge theory (LGT) and the simph- 
cial gauge theory (SGT), we will in this section give a brief overview of the 
discretization procedure from LGT. For a more complete description see e.g. 

m- 

The discretization procedure of both LGT and SGT is based on the following 
identity. Consider a small surface E with area proportional to h^, where h is a 
small positive quantity. Then the following identity holds 

F{A)=n{A)-l + 0{h^), 

s 

where 'H{A) is the holonomy of the one- form A, i.e. the parallel transport 
induced by A around the boundary of E. This parallel transport is defined as 
follows. Given a curve 7 : [0, 1] — >■ M, such that 7(0) — x and 7(1) — y, the 
parallel transport operator along 7 is given by 



U^{x,y)^p(eMij A)] 



where P denotes path-ordering, and the subscript 7 is attached to U to denote 
the path dependence. In LGT, this quantity is known as the Wilson line. 

In LGT, spacetime M is usually discretized by a uniform hypercubic lattice 
L. Neighbouring node positions are related through translation vectors {a^} for 
which we assume \a^\ = h for all /i. To each edge e which connects neighbouring 
nodes, n and n + a^ for some /i, we attach an approximation of the parallel 
transport operator along e. Thus, 

t/^(n) :=exp(i/iA^(n+ -a^)) « C/e(n,n + a^) =Pf exp(i / A) \ . (10) 

In LGT this quantity is called a link variable, link matrix or link group element. 
Furthermore, given a face / of a cube in the mesh, called a plaquette, we 
approximate the holonomy associated to this face as the path-ordered product 
of the link variables along its boundary. In other words, if / lies in the ^v plane, 
with nodes n, n-\- a^, n + a^,, and n + afj_ + a^, we approximate the holonomy as 

Uf{n) := C/^,(n) := U^{n)U,{n + a^)U^(n + a,)U^ [n] 

^ n{A) := P fexp{i (f) A)Y ^^^^ 

where df denotes the boundary of the plaquette /. In LGT, this quantity is 
known as the Wilson loop. Moreover, we approximate the curvature as 

F^u ^Uf-1. (12) 

Finally, the LGT action is defined as 



/ / 



where /3 is related to the coupUng constant by /3 = 4/e^. A discrete gauge 
transformation is associated with a choice of G{n) e SU{2) for each node n. 
Each hnk variable then transforms as 

U^in) ^ G{n)U^in)Gin + a^)-\ (14) 

By the cyclic invariance of the trace, the action Slot is discretely gauge invari- 
ant. 

Remarks. The LGT action can be viewed as a mass lumped FEM action, and 
this observation is useful to have in mind when we construct the simplicial ana- 
logue. In the FEM setting, the gauge potential is assumed to be a lowest order 
curl-conforming Ndlec element in 4d on hypercubes, with one dimension repre- 
senting time [19] . The degree of freedom associated to such a gauge potential 
at an edge e from n to n + a^^ is 






The parallel transport operator is as in equation (10), i.e. U^{n) — exp(iAe) 



Then, the holonomy is approximated as in equation (11), the curvature as in 



equation (12), and one considers [/^ — 1 as the components of the two- form 



where {w/} are the Ndlec basis two- forms. The FEM action associated to such 
a two form is 

2 ^_^, Jm 

where Mff is called the mass matrix, and (•) denotes the scalar product of 
alternating forms w.r.t. the metric. The mass matrix is not diagonal, which 
means that the discrete curvature at different faces interact. This again implies 
that the action is not discretely gauge invariant. However, by diagonalizing the 
mass matrix using numerical quadrature, this action reduces to the LGT action. 



equation ( 13 1. The diagonalization procedure can also be shown to be numerical 



consistent in the sense of approximation theory ^2^ . 

2.3. Simplicial gauge theory 

In this section we construct the discretely gauge invariant simplicial gauge 



theory (SGT) action on a simplicial complex, as defined in appendix Appendix 
[K] The construction is the simplicial analogue of the FEM action described 
above, including additional parallel transport operators to make it discretely 
gauge invariant. 



The curvature associated to the temporal and spatial faces is defined exactly 
as in LGT. In the notation of appendix [Appendix A consider a temporal and 
spatial face 

Mt) := {ir, JT,Jr+At, Jr+At} , , 

/(r) := {ir,jr,kr}, 

where v denotes node i at time t. The time dependency will from here on often 
be suppressed, unless confusion can arise. The spatial and temporal holonomies 
associated to these faces, induced by the gauge potential, are approximated as 

Ufti-ir) = U{ir,JT)U{jr,JT+At)U{jr+At,iT+At)U{ir+At,iT) 
Uf{t) = U{l,j)U{j,k)Uik,l), 

where the arguments v and i are included to indicate where the holonomy is 
located, and the parallel transport operators are defined exactly as in LGT, i.e. 



equation (10). We observe that the holonomies located at different nodes are 
related through the formulas 

Ufti^T + At) = U{ir+At,iT)Uf,{ir)U{ir,iT+At) 

Uf{j)^U{j,{)Uf{i)U{i,j), 

which give formulas for parallel transport of curvature. Hence, we have defined 
the curvature associated to the temporal and spatial faces in our 4d mesh. The 
distinguished point of / and ft, i.e. the location of their holonomy, are denoted 
/ and ft respectively. Note that under a discrete gauge transformation, the 
parallel transport operators are transformed as in LGT, i.e. 

C/(V,V+At) l-> G{ir)U{ir,iT+At)G{ir+At)~^ 

U{i,3)^G{i)U{i,])G{])-\ 

for G{i) e SU{2) for each vertex i. 

As in LGT the curvature is approximated as 

F* « C/r - 1 

(17) 

considered as components of the two-forms 

ft 

J2(Uf^mf, 
f 

where the A are basis functions as described in appendix [Appendix A[ The 
associated FEM action is S — St + Ss, where the temporal part is 

^* = f 3ff 5] M/,/,'tr [{Uf, - l)iUf, - If] , Mf^f, := / A^, • A;-, (18) 

^ , ,, Jm 

It, It 



and the spatial part is 

s^ = fa^E^^'t^ [(^/ - 1)('^/' - ^)"] ' ^ff ■■= [ A/ • V, (19) 

2 ^^^, Jm 

where /? — 2/e^. Note that we have suppressed the dependency of S on A. 
Again, Mf^f and Mff are called mass matrices that depend on the details of 
the mesh, and are described more in detail in appendix |Appendix A[ As pointed 
out in the discussion about the FEM formulation of LGT, the mass matrices 
are not diagonal. This implies that the action is not discretely gauge invariant. 
However, this can be resolved by parallel transport of curvature. The temporal 
and spatial part of the action, St and 5^, are now treated separately. 

The temporal part. Let /((r) and //(t) be two temporal faces. We now use some 
properties of the basis functions, which are explained in appendix [Appendix A 



Since the temporal basis face functions (A/J are piecewise constant in time, 
the interactions between the temporal curvature occur only at coinciding time 
intervals. Also, by properties of the edge basis functions (Ag), which define 
the temporal basis face functions, we can connect the curvature at ft with the 
curvature at // by parallel transport along at most one edge. Thus, we connect 
the curvatures by parallel transport along the connecting edge e = {/t,//} of 
their distinguished points. In other words, we approximate the temporal part 
of the action by 

ft{r)Ji(r) (20) 

X tr [U{j[jt) {Uf^ir) - 1) Uift,n) (C/^.(.) - l)""^ . 

The spatial part. Let / and /' be two spatial faces of a tetrahedron T. The 
curvature associated to the face / at time r will interact with the curvature 
associated to the face /' not only at time r, but also at times r ± At, since the 
facial basis functions are piecewise afSne in time. Thus, to connect the curvature 
at /(t) with the curvature at /'(t') we must parallel transport in both space 
and time. Thus, we replace 

(Ufir) - l)(C//'(.') - 1)^ 

by 

Uifir), /(r))(C//(.)-l)C/(/(r), /'(r))C/(/'(r), f {r'mUf,^^,^-l)"Uif {t')J' (r)) 



in the FEM action (19). In words, we first parallel transport the curvature 
associated to / , located at the vertex /(r)) to the vertex /'(r) along the edge 
e = {/(t), /'(t)}- Then we parallel transport it in the temporal direction from 




Figure 1: Elementary 3d mesh building block containing six tetrahedra, all of 
which share the single interior diagonal. 

/'(r) to /'(r'). So, we approximate the spatial part of the action as 



^SGT 



§3? Y. Mf^r)j'(r')tT [c/(/'(r), /(r)) (C/^(.) - 1) C/(/(r), /'(r)) 



f{r)J'{r') 



H 



xt/(/'(r),/'(r')) (C//-(.o - 1) C/(/'(r'),/'(r)) 



The simplicial gauge theory action is then defined as 

SsGT := SgQj, + Ss(j;rp, 



(21) 



(22) 



and by the cyclic invariance of the trace, this action is discretely gauge invariant. 
A companion paper [22] contains more details about this construction, as well as 
mathematical proofs of consistency with the continuous action ([5]) in the sense 
of approximation theory. 



3. Computer simulation 

For our SGT computer simulations, we chose the euclidean cubic domain 
M = [0, 1]^ C M^ with periodic boundary conditions. We simulated the pure 
gauge SGT action [22] in temporal gauge on a simplicial lattice with the gauge 
group SU{2) . Choice of gauge is not necessary, but it does simplify the algorithm 
slightly, since all temporal edge matrices then reduce to the identity. 

The spatial lattice was constructed using a cubic arrangement of N^ identical 
building block cubes of size h^ , each consisting of six tetrahedra as shown in 
figure [T] The resulting spatial mesh was repeated at TV consecutive time steps 
to form a cubic domain of physical volume (/iTV)^. As described above, each 
spatial edge is part of two temporal square-shaped faces, going forward and 
backward in time. 



The SGT action employs parallel transport matrices in order for gauge in- 
variance to be respected. By defining the distinguished points of all spatial and 
temporal faces to coincide for as many pairs of faces as possible, we only need 
the parallel transport matrices for terms in the action involving pairs of tem- 
poral faces with no common nodes. More details regarding the exact computer 
implementation are given in appendix [Appendix C[ 

3.1. Convergence of the action 

In order to check the continuum limit of the discrete action, we examined 
four different gauge field configurations for which the exact continuum value 
■^cont ^^ ^^^ action is calculable. We did numerical calculations for square 
meshes with A^ = 4, 8, 16, 32 in order to observe convergence of the numerical 
values towards the exact values. By the estimates in [55] we expect that the 
error be of second order in the lattice constant h. We used the following gauge 
field configuration cases: 

1. Gauge field oriented towards the a:-direction in space and towards t^ within 
su(2), with a sinusoidal i-dependence. The only nonzero component of the 
gauge field A is 

A3(t,a;,y,z) := ^sin(27rt), S = I. 

2. Gauge field oriented towards the y-direction in space and t^ within su(2), 
with a sinusoidal a;-dependence. The nonzero component of the gauge field 
in this case was 

Al{t,x,y,z) := —s\n{2-Kx), 5=1. 



3. A case with two nonzero components, 

^i := ^sin(27r2/), A^ := — sin(27rx), S 



2 



27r ' ^" '^ ■ 27r ' " 2 8(27r)4- 

4. A constant field that only contributes to the nonlinear term in the field 
strength, 

-^x '■— Ve, Ay :— Ve, S = -. 

In order to provoke a sizable nonlinear contribution case 3, we chose a small 
/3 — 2/e^ = 1/5. The link matrices needed to evaluate the SGT action are 



calculated from these gauge fields by means of the exponential map (10). 

The results are displayed using double logarithmic plots in figure [2] for tra- 
ditional Wilson action LGT as well as the SGT results. As expected from the 
estimates in p5| , in all cases the relative error behaves as 

Relative error ^ Ch? , 

as determined by extracting the linear coefficient of the second order polynomial 
fits shown in the figures. Note that while the convergence exponent of h is the 



10 




(a) LGT 



Sio- 




(b) SGT 



Figure 2: The relative error of the action versus the number of lattice sites per 



side N, for the actions 1, 2, 3, 4 described in section 3.1 The squares are the 
simulation data points and the solid lines are the second order polynomial fits. 
Errors are proportional to h^ in all cases. 



same in all cases, the prefactor C is smaller in the SGT cases involving time- 
independent fields, due to its finer spatial discretization for the same N. Where 
time-dependence is involved, the errors coincide since the time-discretization 
we have chosen for this SGT simulation is of the same quality as for the LGT 
simulation. 

3.2. Quantum field simulation 

Analogous to the traditional lattice QCD simulations, we performed parallel 
SU{2) quantum field theory Monte Carlo simulations for N = 8. In this case, the 
edge matrices are sampled directly without reference to a gauge field and lattice 
constant value. Therefore, the physical size of the simulation domain is unknown 
prior to experimental comparisons. All dimensional observable quantities are 
automatically calculated in units of powers of the lattice constant h. 

As is customary, it is a Monte Carlo simulation using the Metropolis algo- 
rithm to generate a Markov chain of gauge field configurations that are dis- 
tributed according to the Boltzmann weight exp(— S*). Each Monte Carlo step 
involves randomization of some edge SU{2) matrices, which is done by multi- 
plication of a small su(2) algebra matrix, together with a Metropolis step for 
acceptance/rejection of the update. The algorithm adapted itself to drive the 
MC acceptance rate towards 1/2. Monte Carlo convergence tests were done 
and high quality error estimates were made using data blocking |23j . In addi- 
tion, convergence was verified subjectively by inspection of the time series for 
observable values with their accompanying distributions, as well as time series 
for cumulative averages. The data blocking error estimates were found to be 
smaller than the displayed data points in all the plots. 

We simulated at different values of /3, at each of which we measured the 
average action density S/N^, and a list of different Wilson loops shown in figure 
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Figure 3: The simulated Wilson loops shapes correspond to the outer edges of 
these figures. They lie in the xy^ yz and zx planes. 



[3J all of which are gauge-invariant quantities. For each Wilson loop shape, we 
average over all possible loop positions, as well as loop orientations in the xy, 
yz and zx planes. For a given closed path C, the corresponding Wilson loop 
variable for gauge group SU(n) is defined as 



Wc :-- 



\mv\{u,, 



(23) 



eec 



which involves an ordered product of the edge matrices {Ue\ along the path C. 
Expectation values for any observable quantity O, e.g. the action density 
S/N^ or a Wilson loop Wc: is given by 



{O) 



ly mdC/ejOexp(-5), 



where the partition function Z is defined by 



Z 




exp(-S'). 



(24) 



(25) 



The integration measure involved in these expressions is a product of the nor- 
malized Haar integration measure for each edge group element in the mesh. 
Note that the normalized Haar measure satisfies 



dU = 1. 



(26) 



To accompany these measurements, the strong (small /?) and weak (large 



j3) coupling asymptotic behaviour were calculated in appendix Appendix B 



using methods described in [5]. At strong coupling, this involves various group 
integrals, while at weak coupling it suffices to use a thermodynamic analogy to 
determine the limiting behaviour. 

The simulated results for the action density and Wilson loops are displayed 
in figure |4J In figure |4a| we can see the characteristic and nontrivial behaviour 
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(a) Action density (S'/A''^) versus /3 



(b) Wilson loops {W) versus /3 



Figure 4: P lots showing the /3-dependency of (a) the average action density 
{S/N'^) and (b) the various Wilson loops (W) from figure [sj Solid squares 
are data points and solid lines are linear interpolations. The strong and weak 
coupling asymptotes calculated in [Appendix B| are included for the action 
density and the elementary triangular loop. Monte Carlo errors are smaller 
than the data points. 



in the medium coupling range /? € (1,3). This coincides qualitatively with LGT 
simulations ^^. Only qualitative, not exact, agreement is expected, since the 
physical lattice constant will differ in each type of simulation. Compared to 
LGT simulations, the behaviour at small /3 deviates more from linearity due 
to the nonlinear aspects of the SGT action. In this region, the actions do not 
approximate the continuum action, and differences between discrete actions are 
unphysical. 



The Wilson loops in figure 4b show the same qualitative behaviour as do 
LGT simulation results, and approaches the calculated asymptotes nicely. Also 
here, the behaviour is less linear at small /3 for the same reason as stated above. 
The typical strong suppression of the Wilson loops as functions of loop area is 
reproduced, as expected from the area law behaviour that indicates confinement. 



4. Conclusions 

We have implemented the general SGT action on a particular simplicial 
mesh, and performed Monte Carlo quantum field theory simulations that show 
sensible results that are qualitatively consistent with standard LGT simulations, 
as must be the case for this initial proof-of-concept implementation. 

We expect that this method will lend itself nicely to the use of mesh refine- 
ment within quantum QGD simulations, and that this will lead to opportunities 
of novel applications using nontrivial mesh structures, e.g. in the vicinity of 
gluon flux tubes as mentioned in the introduction. 

The nondiagonal nature of the action increases the amount of computer work 
in the Metropolis step after each proposed update. However, since the number 
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of interactions for each elementary face is finite, the scahng at large meshes for 
this model will be the same as for traditional QCD. There might be possibilities 
of real-time adaptive diagonalization, thereby increasing the algorithm efficiency 
throughout the initial part of the simulation. 

Appendix A. Simplicial complex, finite elements and mass matrices 

Consider a collection of vertexes, edges, faces, tetrahedra in 3d space. These 
elementary objects are called simplexes, and the collection of these a simplicial 
complex T. For any /c-dimensional simplex Tk for 1 < fc < 3, the boundary 
dTk is a union of (k — l)-dimensional simplexes. Consult [25' Section 5.1] for 
a precise definition. In our construction, we assume that this spatial simplicial 
complex spans the spatial domain S. The vertexes, edges, faces, and tetrahedra 
according to dimension, and are labeled i, e, /, and T respectively. The symbol 
T will be used for simplexes of any dimension. 

In order to expand this to a 4d spacetime simplicial complex T, consider a 
uniform time-discretization with a time-spacing At. The simplicial complex T 
is then repeated at each discrete time step value r. For each such r, we define 
additional simplexes for our T by extruding each simplex of T along the time 
interval [r, r -I- At] . As the basic building block in classical 3d FEM theory 
is a tetrahedron T, the basic building block in this extended FEM version is 
T X /t, where I-r = [t,t + At], i.e. a time-extrution of a tetrahedron. Temporal 
edges are generated by extruding 3d vertices, and temporal faces by extruding 
3d edges. 

The space of Whitney k-forms on T (T) is denoted W^{T) {W^{T)), with 
canonical basis (At), T ranging over the set of fc-dimensional simplexes in 7" 
[20j . The 0-forms Xi are the barycentric coordinate maps for each vertex i. In 
other words, it is the piecewise affine map taking the value 1 at the vertex i 
and at other vertices. For an edge e = {i,j}, with orientation i -^ j, the 
associated Whitney 1-form is defined by 

Ae := Ay := A.VA^- - A^ VA,. (A.l) 

For a face / = {i,j, k}, whose orientation is i — > j — > fc, the associated Whitney 
2-form is defined by 

Xf := A,yfe := 2 (A.VA^ x VAj, + Aj VA^ x VA, + AfeVA, x VA^) . (A.2) 

In the 4d spacetime FEM setting, these basis fc-forms are extended to be 
piecewise affine in time and are denoted {Ax(t)), i-C 

At -^ t^T(T) = Xt<^PI, 

where P* denotes polynomials in the time variable of degree at most one, and 
T{t) := (r, T) denotes the spatial simplex T at temporal node r. More precisely, 
At(t) is the piecewise affine function in time, taking the value At at t and at 
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the other temporal nodes. In addition, we define temporal basis edge and face 
functions. 

To every vertex i in the spatial mesh there are temporal edges etir) = 
{VjV+At}; where v :— ^(t). The temporal basis edge function attached to 
et(T) is then the piecewise constant function in time defined by 

iK..i,dt, te[r,r + At] 
10, otherwise. 

where tt is the canonical projection onto the space S, 

t: ■.M = Rx S ^ S, 

and dt is the standard basis one-form in the temporal direction. 

To every spatial edge e there are corresponding temporal faces /((t) = ex/,-. 
The temporal basis face function attached to ft (t) is then the piecewise constant 
function in time defined by 

(Ko.A^^dt, te[r,r + At] 
^'^ '^ ' [0, otherwise. 

In addition to these basis functions, we must define mass matrix elements. 
Let rriTT' denote the classical 3d mass matrices for spatial Whitney elements 

TflTT' -- I ^T ' ^T' ; 

Js 

where T, T' are fc-dimensional simplexes, and (•) denotes the scalar product of 
alternating forms. 

In the definition of the SGT action we use the generalization 

MT(t)T'{T) ^ / AT(t) • At'(t)- 

This generalization can be expressed through the classical mass matrices by 
performing the time integration explicitly. Thus, let T be a spatial tetrahedron 
and Ir = [t,t + Ai] . Considering now only this time interval, the piecewise 
affine function taking the value 1 at time t and at time r + Af is given by 

, , i — T 

Pr{t) = l-^. 

The analogous function for the temporal node r + At on the same time interval 
is given by 

Restricted to the basic building block T x /t-, we therefore get 
^/(r)/'(r)(T X /O = / A/(,) • A/,(,) ^ f pI f Xf Xf, ^ ^Atmff,iT), 

JTxI-r Jl^ JT -J 
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Mf(r)f'iT+At)iT X Ir) - / Ay(^) • Aj,(^+At) 

JTxI^ 



Mf 



= / PrPr+At / A/ • A// = -Atmff,{T), 
(r+At)f'{T+At){T X /t) = / Ay(T._|_At) ■ A/'(r+At) = 

Similarly, the mass matrix element corresponding to the temporal face basis is 
given by 

MMr)fiir){T Xlr)^ f A/,(.) • Af^f^r) ^ ^J ^e ' K ^ 7^"^ee'(T). 



Appendix B. Strong and weak coupling limits 

Appendix B.l. Strong coupling limit 

Here we will show some details regarding the calculation of the strong cou- 
pling limits of the elementary triangular Wilson loop. We will use the following 
integrals over SU{2) group space [Sj 

f dUU"'^ = 0, f dUU^^f^^U""^^"^ = ijai«2^ft&^ 
■' ■' . \ (B.l) 

where the Greek symbols are matrix indices. 

In this calculation, the Wilson loop encircles an elementary spatial triangular 



plaquette Pt at time t. We denote this Wilson loop by Wp^. By equation (23), 
it is given by 

Wp^ := ^mr (UaUbUc) , 

where the plaquette Pt is encircled cyclically by the SU{2) edge matrices Ua,Ub 
and Uc- Due to our choice of distinguished points and plaquette orientations, 
the spatial SGT action is given by 

^ fj' 
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where the sum extends over ah spatial faces at all times. Since we are interested 
in small /3, consider a first order truncated Taylor expansion of the exponential 



in equation (24), i.c 




mr (UaUbUc) Y, Mff'tT {UfUf, ~Uf- Uf, 



fj' 

By the properties of the SU(2) integration measure, terms involving integration 
over odd powers of link matrices vanish. Therefore, nonvanishing contributions 
to the integral only come from terms where either / and/or /' coincide with 
the plaquette Pt. The UfU¥, doesn't contribute. Indeed, if either / or /' differ 
from Pt, we such a term includes an integral over a single power, which vanishes. 
If on the other hand / = /' = Pt, we have UfU¥ = 1 which again leads to 
an integral over a single power and thus vanishes. This is also the case for the 
constant term in the parenthesis. 
We are left with 

(Wp,) « ^n I [\{duA tr (C/,[/,[/e) Y. ^^ff'^' (Uf + Up) , 

1^ -' \ e / fJ' 

where we have moved the real part operator 5ft outside of the integral. Contri- 
butions only come when at least one of /, /' coincide with Pt. Therefore, by the 
properties of the particular mesh we have constructed, 

(Wp,) - J- {Mp^P, + Mp^p^^, + Mp,P,_,) X 



x5ft 



J I Y[dU A tr (UaUhU,) tr {Up^ + C/^) 



Using C/pj := UaUbUc and the SU{2) integration formulas (B.l), we get 



{Wp,) « ^ {Mp,p, + Mp,p,^, + Mp,P,_,) = p, 

where we have used Z « 1 for small /3. The last equality follows from the par- 
ticular mass matrix element values produced by our choice of simplicial lattice. 
A similar calculation, only slightly more involved because several faces are 
involved, can be performed to determine the strong coupling limit of the action. 
Approximations of higher order in /3 can be found by including higher order 
terms in the Taylor expansion of the exponential. 

Appendix B.2. Weak coupling 

In order to determine the weak coupling limit of the action density, we simple 
follow a thermodynamic analogy described in [3]. At large /3, the system is 
described well by a gaussian partition function approximation. This corresponds 
to a free theory, and we can find the weak coupling limit of the action by 
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distributing an amount kT/2 = 1/2/9 of energy among all the degrees of freedom 
in the theory. We have seven edges for each building block cube, each of which 
contributes three degrees of freedom (the number of generators of SU{2)). To 
obtain the action, we multiply by /3, which results in 

SsGT ^ /3 X ^ X 7 X 3 = y 7V^ as /3 ^ cx). (B.2) 

This result can be used to determine the same limit of the triangular Wilson 
loop in the Q;/3-plane. We have 

where there is no sum over the spacetime indices. The antisymmetric field 
strength has six independent spacetime components. By the equipartitioning of 
the euclidean energy among these degrees of freedom, we have 

/frrp2 \ _ l/.,/p pM^u _ V/ -^SGT \ _ 42^2 



Now using /3 = 2/g^, we get 



(^i) = l-^- (B-3) 



Appendix C. Computer implementation 

Our computer implementation of the simplicial lattice and accompanying 
SGT action consists of object-oriented C+-I- code, using MPICH2 [35] for par- 
allelization, running on a quadruple CPU run-of-the-mill modern workstation 
computer. The data structures involved are reminiscent of what is used in im- 
plementations of the finite element method. This involves different types of 
mass matrix and connectivity information for elements of the simplicial mesh. 
The parallelization consisted of running independent simulations on each node, 
and averaging the results. We used the yarn2 algorithm from the TINA pseudo- 
random number generator |27) , which is designed for use in parallelized algo- 
rithms. Although the edge matrix randomization appeared to perform stably 
enough for our purposes, we regularly did projections of the edge matrices onto 
SU{2) as a precautionary measure. 
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